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Motivated by complex multi-fluid geometries currently being explored in fibre-device 
manufacturing, we study capillary instabilities in concentric cylindrical flows of N fluids 
with arbitrary viscosities, thicknesses, densities, and surface tensions in both the Stokes 
regime and for the full Navier-Stokes problem. Generalizing previous work by Tomotika 
(N = 2), Stone & Brenner (N = 3, equal viscosities) and others, we present a full 
linear stability analysis of the growth modes and rates, reducing the system to a linear 
generalized eigcnproblcm in the Stokes case. Furthermore, we demonstrate by Plateau- 
style geometrical arguments that only axisymmctric instabilities need be considered. We 
show that the N = 3 case is already sufficient to obtain several interesting phenomena: 
limiting cases of thin shells or low shell viscosity that reduce to N = 2 problems, and 
a system with competing breakup processes at very different length scales. The latter 
is demonstrated with full 3-dimensional Stokes-flow simulations. Many N > 3 cases 
remain to be explored, and as a first step we discuss two illustrative N — > oo cases, an 
alternating-layer structure and a geometry with a continuously varying viscosity. 
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1. Introduction 

In this paper, we generalize previous linear stability analyses ( |Lord Rayleigh 1879 
|1892| |Tomotika||T935l |Stone fc Brenner||1996| |Chauhan erZ|2000| ) of Plateau-Rayleigh 
(capillary) instabilities in fluid cylinders to handle any number (N) of concentric cylindri- 
cal fluid shells with arbitrary thicknesses, viscosities, densities, and surface tensions. This 
analysis is motivated by the fact that experimental work is currently studying increas- 
ingly complicated fluid systems for device-fabrication applications, such as drawing of 



microstructured optical fibres with concentric shells of different glasses/polymers (Hart 


et a/.|2002 


Kuriki et a/.|2004 Pone et al. 2006; Abouraddy et aL|2007||Sorin et aZ.|2007 


Deng et al. 


2008 1 or generating double emulsions ( 


Utada et al. 


2005 Shah et al. 


2008). 



Although real experimental geometries may not be exactly concentric, we show that sur- 
face tension alone, in the absence of other forces, will tend to eliminate small deviations 
from concentricity. We show that our solution reduces to known results in several limiting 
cases, and we also validate it with full 3-dimensional Stokes-flow simulations. In addition, 
we show results for a number of situations that have not been previously studied. For the 
limiting case of a thin shell, we show a connection to the classic single-cylinder and flat- 



plane results, consistent with a similar result for air-clad two-fluid jets ( Chauhan et al. 
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2000). In a three-fluid system, we exhibit an interesting case in which two growth modes at 
different wavelengths have the same effective growth rate, leading to competing breakup 
processes that we demonstrate with full 3-dimensional Stokes-flow simulations. We also 
consider some many-layer cases, including a limiting situation of a continuously varying 
viscosity. Using a simple geometrical argument, we generalize previous results ( |Plateau| 
1873 Lord Rayleigh|1879 Chandrasekhar|1961 ) to show that only axial (not azimuthal) 



instabilities need be considered for cylindrical shells. Numerically, we show that the sta- 
bility analysis in the Stokes regime can be reduced to a generalized eigenproblem whose 
solutions are the growth modes, which is easily tractable even for large numbers of layers. 
Like several previous authors (Tomotika 1935 |Stone fc Brenner 1996 Gunawan et al. 



2002 2004), we begin by considering the Stokes (low- Reynolds) regime, which is con- 



2007 



sistent with the high viscosities found in drawn- fibre devices ( Abouraddy et al. 
Deng et aZ.||2008 1. In Appendix [C] we generalize the analysis to the full incompressible 
Navier-Stokes problem, which turns out to be a relatively minor modification once the 
Stokes problem is understood, although it has the complication of yielding an unavoid- 
ably nonlinear eigenproblem for the growth modes. Semi-analytical methods are a crucial 
complement to large-scale Navier-Stokes simulations (or experiments) in studying capil- 
lary instabilities, since the former allow rapid exploration of wide parameter regimes (e.g. 
for materials design) as well as rigorous asymptotic results, while the latter can capture 
the culmination of the breakup process as it grows beyond the linear regime. 



Capillary instability of liquid jets has been widely studied (|Lm||2003 Eggers & Villcr- 
since 



maux 



2008) 



Plateau (1873) first showed that whenever a cylindrical jet's length 
exceeds its circumference, it is always unstable due to capillary forces (surface tension). 
Plateau used simple geometrical arguments based on comparing surface energies before 
and after small perturbations. Lord Rayleigh introduced the powerful tool of linear sta- 



bility analysis and reconsidered inviscid water jets ( |Lord Rayleigh 1879 ) and viscous 
liquid jets (Lord Rayleigh 1892). In linear stability analysis, one expands the radius R 

R + 8Re ikz ~' lut , where 5R < R , 
is an exponential growth rate. 



as a function of axial coordinate z in the form R(z) 
2n/k is a wavelength of the instability, and a — Imw 
Given a geometry, one solves for the dispersion relation(s) w(fc) and considers the most 
unstable growth mode with the growth rate tr max to determine the time scale of the 
breakup process. The wavelength 27r/A; corresponding to cr ma x has been experimentally 



verified to match the disintegration size of liquid jets (Eggers & Villermaux 20081. By 
considering the effect of the surrounding fluid, Tomotika ( 19351) generalized this analysis 



to a cylindrical viscous liquid surrounded by another viscous fluid, obtaining a 4 x 4 
determinant equation for the dispersion relation. (Rayleigh's results are obtained in the 
limit of vanishing outer viscosity.) A few more limiting solutions have been obtained 
(Meister & Scheclc 1967 Kinoshita et a/.|[T994[ ) by generalizing Tomotika's approach to 
non-Stokes regimes. Beyond the regime of a single cylindrical jet, Stone & Brenner ( 1996 ) 
analysed the three-fluid (N = 3) Stokes cylinder problem, but only for equal viscosities. 
Chauhan et al. ( 2000 ) analysed the N = 3 case where the inner two fluids have arbitrary 



viscosities and the outermost fluid is inviscid gas, taking into account the full Navier- 



Stokes equations. Gunawan et al. (2002 2004) considered an array of identical viscous 



cylinders (in a single row or in a triangular configuration). 

Much more complicated, multi-fluid geometries are now being considered in experi- 



mental fabrication of various devices by a fibre-drawing process ( 


Hart et al. 2002 


Kuriki 


et aZ.||2004 


Pone et aZ.|2006||Abouraddy et aZ.||2007[ Sorin et a/.|2007 


Deng et al. 


2008). 



In fibre drawing, a scale model (preform) of the desired device is heated to a viscous state 
and then pulled (drawn) to yield a long fibre with (ordinarily) identical cross-section but 
much smaller diameter. For example, concentric layers of different polymers and glasses 
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can be drawn into a long fibre with submicron-scale layers that act as optical devices for 
wavelengths on the same scale as the layer thicknesses (Abouraddy et al. 2007). Other 
devices, such as photodetectors ( Sorin et a/.|2007 ), semiconductor filaments (Deng et al. 



2008 2010), and piezoelectric pressure sensors (|Egusa et al. 2010) have similarly been 



incorporated into microstructured fibre devices. That work motivates greater theoretical 
investigation of multi-fluid geometries, and in particular the stability (or instability time 
scale) of different geometries is critical in order to predict whether they can be fabri- 
cated successfully. For example, an interesting azimuthal breakup process was observed 
experimentally ( Deng et aZ.||2008 ) and has yet to be explained (?); however, we show in 
this paper that azimuthal instability does not arise in purely cylindrical geometries and 
must stem from the rapid taper of the fibre radius from centimetres to millimetres (the 
drawn-down "neck"), or some other physical influence. [For example, there may be clastic 
effects (since fibres are drawn under tension), thermal gradients at longer length scales 
in the fibre (although breakup occurs at small scales where temperatures are nearly uni- 
form), and long-range (e.g., van der Waals) interactions in submicron-scale films.] Aside 
from fibre drawing, recent authors have investigated multi-fluid microcapillary devices, 
in which the instabilities are exploited to generate double emulsions, i.e. droplets within 
droplets ( Utada et q/.|2005 ). Because the available theory was limited to equal viscosities, 
the experimental researchers chose only fluids in that regime, whereas our paper opens 
the possibilities of predictions for unequal viscosity and more than three fluids. 

We now formulate the mathematical problem that we solve, as depicted schematically 
in figure [I] The total number of viscous fluids is N and the viscosity of the n-th (n = 
1,2, ••• ,N) fluid is ^W. The surface-tension coefficient between the n-th and (n+1) — th 
fluid is denoted by 7W. For the unperturbed steady state (figure lb), we assume that 
the n-th fluid is in a cylindrical shell geometry with outer radius and inner radius 
R(n-i) < The first (n = 1) fluid is the innermost core and the N-th fluid is the 

outermost one (extending to infinity), so we set i? (0) = and i? (JV ) = +00. To begin 
with, this system is analysed in the Stokes regime (low Reynolds number) and we also 



neglect gravity [in the large Froude number limit, valid for fibre-drawing (Deng et al. 
2011)], so the fluid densities are irrelevant. In Appendix [Cj we extend this analysis to 



the Navier-Stokes regime, including an inertia term for each layer (with density p^ 71 '). 
As noted above, linear stability analysis consists of perturbing each interface R^ by a 
small sinusoidal amount §R( n ) e lkz - luJt ( to lowest order in 5R( n \ Stokes' equations are then 
solved in each layer in terms of Bessel functions, and matching boundary conditions yields 
a set of equations relating uj and k. Although these equations can be cast in the form of 
a polynomial root-finding problem, similar to Tomotika (1935), such a formulation turns 
out to be ill-conditioned for large N, and instead we formulate it in the Stokes regime 
as a generalized eigenproblem of the form M2(fc)£ = iu;Mi(fc)£, which is easily solved 
for the dispersion relations w(fc) (with the corresponding eigenvectors £ yielding the 
relative amplitudes of each layer). In the Navier-Stokes regime, this becomes a nonlinear 
eigenproblem. 



2. Azimuthal stability 

For any coupled-fluids system of the type described in figure [T] a natural question to 
ask is whether that system is stable subject to a small perturbation. If an interface with 
area S has surface energy jS, then a simple way to check stability is to compare surface 
energies (areas) for an initial configuration and a slightly deformed configuration with 
the same volume. In this way, it was shown that any azimuthal deformation is stable for 
a single cylindrical jet (ILord Rayleigh||1879| |Chandrasekhar||1961|). Here, we employ a 
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(a) Azimuthal cross-section (b) Steady state (c) Perturbed state 



Figure 1: Schematic of the concentric-cylinder geometry considered in this paper, (a) 
Cross-section of N layers and corresponding radii R^ n \ viscosities and surface ten- 
sions y( n K Starting with the perfect cylindrical geometry (b), we then introduce small 
sinusoidal perturbations (c) and analyse their growth with linear stability analysis. 



similar approach to demonstrate that the same property also holds for multiple concentric 
cylindrical shells. Note that this analysis only indicates whether a system is stable; in 
order to determine the time scale of an instability, we must use linear stability analysis 
as described in subsequent sections. 

For the unperturbed system, we define the level-set function <p") = r — R^ n ' . < -"- ) = 
defines the interface between the n-th and (n + l)-th fluids. Similarly, we define the level- 
set function for the perturbed interface (figure lcl between the n-th and (n + l)-th fluids 

by 

0<*>(r,z,(« = r-C (n W). (2.1) 



Following the method of normal modes (Drazin & Reid 2004), in the limit of small 
perturbations, a disturbed interface £^ can be chosen in the form 

C [n] {z,4>) = R {n) +5R {n) e i{kz+m ^ +0[{5R {n) ) 2 ], (2.2) 

Assuming incompressible fluids in each layer so that volume is conserved (and assuming 
that the cylinder is much longer than its diameter so that any inflow/outflow at the end 
facets is negligible), we obtain a relation between f W and R^ 



4i?(") 



■ 0[{5R {n) f] 



(2.3) 



Let S(<j)( n ' ) de note the s urface area of (p n '{r,z,<j)) = in one wavelength 2n/k. From 
equations (2.1) and (2.2), S((f)( n ') can be expressed in cylindrical coordinates as: 



/ <9C (n) 1 9C (n) 

l + (4-) 2 + ( 77 ^4^) 2 d0dz. 



dz ' ' v C (n) 9, 

Now we can compare the total interfacial energy between the unperturbed system £ 



(2.4) 
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(4>W) and the perturbed system £ = Y,n=i 7 (n) S(<£ (n) ): 

N-l 

71=1 



(2.5) 



From the surface-energy point of view, small perturbations will grow only if £ — £ < 
0. Therefore, from (2.5), we can conclude that all the non-axisymmetric perturbations 



(to 7^ 0) will be stable. There is one special case that needs additional consideration: if 
k = and m = 1 in (2.51, the first term is zero, so one must consider the next-order 
term in order to show that this case is indeed stable (i.e. elliptical perturbations decay). 
Even more straightforwardly, however, k — corresponds to a two-dimensional problem, 
in which case it is well known that the minimal surface enclosing a given area is a circle. 



3. Linear stability analysis 

In the previous section, we showed that only axisymmetric perturbations can lead to 
instability of concentric cylinders. Now we will use linear stability analysis to find out 
how fast the axisymmetric perturbations grow and estimate the break up time scale for 
a coupled N-l&yer system. 

Here, we consider fluids in the low-Reynolds-number regime [valid for fibre-drawing 



( Abouraddy et al. 2007 Deng et al. 


2008 


|] and thus the governing 


each fluid are the Stokes equations 


Ockendon & Ockendon|1995 ) 



of the n-th fluid is — [u r n \r, z, i), iti™''^, z,t)}, where ui^ is the radial component 
of the velocity and u^ is the axial component of the v elocity. The dyna mic pressure in 
the n-th fluid is denoted by p( n K The Stokes equations ( Batchelor||l973 ) are 



,0) 



d 2 ul n) 


1 dui n) 


(«) 

Ur 


d 2 U r n) \ 




dr 2 


r dr 




f dz 2 J 


dr 




'd 2 ui n) + 


i dui n) 


d 2 ui n) \ 


dpi") 


dr 2 


r dr 


f dz 2 j 


dz 



and the continuity equation (for incompressible fluids) is 



du 



,(») 



du y z 



O) 



dr 



dz 



0. 



(3.1a) 
(3.16) 

(3.2) 



3.1. Steady state 

Because of the no-slip boundary conditions of viscous fluids, without loss of generality, 
we can take the equilibrium state of the outermost fluid to be 



-,(N) _ 







u 



(JV) _ 



= 



P 



(JV) _ 



= 



(3.3) 



for r > : \ and of the n-th (n < N) fluid to be 
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for R( n -^ < r < (In Appendix [cj we generalize this to nonzero initial relative ve- 

locities for the case of inviscid fluids, where no-slip boundary conditions are not applied.) 



3.2. Perturbed state 

The perturbed interface, corresponding to the level set = r — £ 
axisymmetric perturbation, takes the form 

( {n) (z,t) = R<ri +6R^e i{kz -^ + 0[(<Si? (n) ) 2 ] 

Similarly, the perturbed velocity and pressure are of the form 



(n) 



0, with an 



(3.5) 



ui n \r, z,t) 




"-(«) 

Ur 




~6ul n) (rj 


(™) 1 4-\ 

Uz (r,z,t) 




_(n) 
U z 


+ 


5u { z n \r) 


p(")(r,z,t) 




p(n) 




Sp^{r) 



Note that the Stokes equations (3.1) and continuity equation (3.2) imply that 

Q2 p (n) 



dr 2 



1 dp^ d 2 p^ 
r dr dz 2 



0. 



(3.6) 



(3.7) 



Substituting the third row of (3.6) into (3.7), we obtain an ordinary differential equation 
for SpW (r) 

^ + ~-tAsp^(r) = 0. (3.8) 
v ar z r dr J 

Clearly, 8p^{r) satisfies the modified Bessel equation of order zero in terms of kr. 
Therefore, we have 



( 5p(")(r) = 4 n) X (fcr) + 4" ) 



Io(kr), 



(3.9) 



where Kq(-) and Iq(-) are modified Bessel functions of the first and second kind (Kq is 
exponentially decreasing and singular at origin; I Q is exponential ly g rowing), and c^ n ' 



and are constants to be determined. Inserting 5p^(r) into (3.1) and solving two 



inhomogeneous differential equations, we obtain the radial component of velocity 



M _ -(») rK (kr) {n) rl (kr 

6U r W" C 1 2m („) + C 2 2m („) 



) ( n )Ki(kr) ( n ) h(kr) 



2/j,( n )k 



(3.10) 



and the axial component of velocity 
6u<?\r) =4 n) 



iK (kr) 
jjL^k 



irKi(kr) 
2/x(") 



„(") 



Ho(kr) 
jjL^k 



irli{kr) 



(„)iif (fcO , ( n )i/o(fc^) 



2/i(»)jfe 



(3.11) 



where 



» 



and 



» 



are constants to be determined. Imposing the conditions that the 



velocity and pressure must be finite at r = and r = +oo, we immediately have 



0. 



(3.12) 



3.3. Boundary conditions 

In order to determine the unknown constants c^"' = (c^, <4 , c|j , ) in each layer, 
we close the system by imposing boundary conditions at each interface. Let be the 
unit outward normal vector of interface r = ^ n \z, t) an d t*-™- 1 be t he unit tangent vector. 
Formulae for and i*-™* 1 are given by equations (All and (A2) of Appendix [Xj First, 
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the normal component of the velocity is continuous at the interface, since there is no 
mass transfer across the interface, and so 



u 



(n)| 



\r=C 



(») = u ("+i) 



(") 



(3.13) 



For the at-rest steady state (3.3) and (3.4), this condition is equivalent (to first order) 
to the continuity of radial velocity: 

u^{ r> z,t)\ r=(M =4 n+1 \r,z,t)\ r=cM . (3.14) 

Second, the no-slip boundary condition implies that the tangential component of the 
velocity is continuous at the interface: 



u (n) | r=C (n> -t (n) 



,("■+1)1 



=C<») 



t 



(n) 



(3.15) 



(The generalization to inviscid fluids, where no-slip boundary conditions are not applied, 
is considered in AppendixjCj) For the at-rest steady state ( 3.3 ) and ( 3.4 ), this is equivalent 
(to first order) to the continuity of axial velocity: 

u^(r,z,t)\ r=CM =u^ +1 \r,z,t)\ r=c , ny . (3.16) 

Third, the tangential stress of the fluid is continuous at the interface. The stress tensor 
of the n-th fluid in cylindrical coordinates for axisymmetric flow ( Kundu &: Cohen|2007 1 
can be expressed as 



-O) 



» 



(n) 



dz 



- p (n) 



2[i 



(") 



dr 

~d~T 



dz dr 

The continuity of the tangential stress at the interface implies that 
n 



rW| r=c( »)-tW 



The leading term of (|3.18| leads to 

'du r n) 



(n) 



du 



(«)' 



dz 



dr 



r=(W — M 



<9< +1 <9< +1 

— 1 — 

dz dr 



\r=((™)- 



(3.17) 



(3.18) 



(3.19) 



Fourth, the jump of the normal stress across the interface must be balanced by the 
surface-tension force per unit area. The equation for normal stress balance at the interface 
is 



(n) . ( T (n+1) _ r (n)\ 



| r=f( „) • n 



(*») — j(n) K (n) 



(3.20) 

where (r, z, t) is the mean curvature of the interface. The curvature can be calculated 
directly from the unit outward normal vector of the interface by = V • (see 
Appendix |A|. Substituting (3.17), (Al), and (A5) into (3.20), we have the following 
equation (accurate to first order in 5R (n ^) for the dynamic boundary condition: 

t du r n) ' 



dr 



= 7 W 



-P 



(«) 



2[i 



(n)'_ 



1 



dr 

6u r n) (R^) 



k A - 



1 



(i?(«)) 2 



(3.21) 



3.4. Dispersion relation 



Subs tituti ng ( |3i3| and (|379[)-( |3. 1 1 1 ) into the boundary conditions [( |3.14[ ), ( |3.16[ ), ( |3.19 l 
and (3.21)] and keeping the leading terms, we obtain a linear system in terms of the 
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unknown constants = (4™ , 4 , 4 » 4 "'). After some algebraic manipulation, these 
equations can be put into a matrix form: 



[ fi(n,n) _|_ g(n) j c (n) _ ^(n,n+l) c (n+l) _ g 



(3.22) 



where /*("•>»') and B( n ) are 4x4 matrices given below. 



A" 



(n) 



2A<" ) -fcfl'")A-<" ) 2J^ ) +fc J R<">J<" ) 
kR (n) K (n) _ R (n) kR (n) jM + j(r 



it"') 



A', 



(n) 



O) 



g(n) = 



7 (™)/fc(l 



(fefl(">) 









(n) 



'1 



(n) 









kR^K^ kR^I { n) K[ r 



rin) 

'o 







r(«) 



,(«> 
kR("> 



(3.23) 
(3.24) 



Here, -Kg , if}™ , and denote the corresponding modified Bessel functions eval- 
uated at kR( n \ Combining the boundary conditions from all N — 1 interfaces, we have 
the matrix equation 

M iN) £ = (3.25) 



for the undetermined constants £ = (4 > , , 

1 -iu> 2 

',4(1,1) _^(l,2) 
/\(2,2) 



^-i), 4^,4^), With 



-^(2,3) 



1 

— iw 



£(N-2,N-2) _^{N-2,N~1) 

MN-X,N-1) _A(N-1,N) 

bw 

e( 2 ) 

e(w-2) 

e(jv-D o 



, (3.26) 



where /U 1 ' 1 ) and are the second and fourth columns of Z^ 1,1 ) and f?W, and /U^ -1 '^ 
i s the first and third columns of A^ N ~ 1,N K To obtain a non-trivial solution of equation 
(3.25), the coefficient matrix M^' must be singular, namely 



det(MW) = 0. 



(3.27) 



Since is zero except in its fourth row, w only occurs in the 4 th, 8th, • • • , 4(N — l)-th 
rows of M( N K The Leibniz formula implies that equation (3.27) is a polynomial in 1/uj 



with degree N— 1. Therefore, we could obtain the dispersion relation lo = u>(k) by solving 



the polynomial equation (3.27). Instead, to counteract roundoff-error problems, we solve 



the corresponding generalized eigenvalue problem as described in §|4] 
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3.5. Eig en- amplitude and maximum growth rate 
The N - 1 roots of (13.27} are denoted by Uj(k), where j = 1,2, • • • , JV - 1. Since M'"' 



is singular when w 

constant. Therefore, for each ujj(k), the corresponding perturbed amplitude SR^ 1 ' on 
the n -th int erfac e can also be determined up to a proportionality constant by equations 



£ can be determined by equation (|3.25|) up to a proportionality 

(n) 



(2) 



(n) 



(N- 



( |3.10p and (|A4j). Let us call the vector 8R 3 = (S Rf\ S R,f } , ■ ■ ■ ,SR^\--' ,SR^ N 
where we normalize = 1, the "eigen-amplitudes" of frequency u)j. Any arbitrary 

initial perturbation amplitudes A = {A^ X \A^ 2 \ ■ ■ ■ ,A( N ~ 1%> ) can be decomposed into a 
linear combination of eigen-amplitudes, namely A = ^2j = i o,j5Rj for some aj. Since 
the whole coupled system is linear, the small initial perturbation Ae lkz will grow as 



The growth rate for a single mode is o~j(k) = Imwj(fc) since our time dependence is 
e -iujt_ jf a . > o, then the mode is unstable. As described above, for an TV-layer system, 
we have N — l different growth rates for a single k, and we denote the largest growth rate 
by Cmax(fc) = max, [aj (k)] . Moreover, the maximum growth rate over all k is denoted by 
°max = m a x fc[c m ax(&)] = max j,fc [<Tj (k)] , and the corresponding eigen-amplitudes of this 
most-unstable mode arc denoted by <5i? max ; fc max denotes the corresponding wavenumber. 



4. Generalized eigenvalue problem 

It is well known that finding the roots of a polynomial via its coefficients is badly 
ill-conditioned ( [Trefethen fc Bau|1997 1 . Correspondingly, we find that solving the deter- 
minant equation (3.27) directly, by treating it as a polynomial, is highly susceptible to 
roundoff errors when N is no t small. In particular, it is tempting to use the block struc- 
ture of to reduce (3.27) to a 4 x 4 determinant problem via a recurrence. However, 
the entries of this 4x4 matrix are high-degree polynomials in uj whose coefficients thereby 
introduce roundoff ill-conditioning. Instead, we can turn this ill-conditioned root-finding 
problem into a generalized eigenvalue problem by exploiting the matrix-pencil structure 
of M v : 



AfW£ = o- M^ N) {k)£ = iojM[ N) (k)£ 



(4.1) 



Thus, finding the dispersion relation uj(k) turns out to be a generalized eigenvalue prob- 



lem with matrices (Mi , : ). Since MY"' is non-singular, this (regular) generalized 
eigenvalue problem is typically well- conditioned ( Demmel fc Kagstromj 1987[ ) and can 
be solved via available numerical methods (Anderson et al. 1999). In principle, further 
efficiency gains could be obtained by exploiting the sparsity of this pencil, but dense 
solvers are easily fast enough for N up to hundreds. 



AN) 



5. Validation of our formulation 

As a validation check, our iV-layer results can be checked against known analytical 
results in various special cases. We can also compare to previous finite-element calcula- 
tions peng et «Z.||2011| >. 



5.1. Tomotika's case: N = 2 



Tomotika| ( |1935 ) discussed the instability of one viscous cylindrical thread surrounded 



by another viscous fluid, which is equivalent to our model with N = 2. It is easy to verify 
that det(M^) = 0, where M< 2 > = [/4« + ^B^, D«], gives the same determinant 



equation as (34) in Tomotika (1935). 
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In contrast with the Stokes-equations approach, Tomotika began with the full Navier- 
Stokes equations, treating the densities of the inner fluid p' and the outer fluid p as small 
parameters and taking a limit to reach the Stokes regime. However, special procedures 
must be taken in order to obtain a meaningful determinant equation in this limit, be- 
cause substituting p' — and p = directly would result in dependent columns in the 
determinant. Tomotika proposed a procedure of expanding various functions in ascending 
powers of p and p', subtracting the leading terms in dependent columns, dividing a quan- 
tity proportional to pp' , and finally taking a limit of p — > and p' — > 0. We generalize 
this idea to the iV-shell problem in Appendix [Cj but such procedures are unnecessary if 
the Stokes equations are used from the beginning. 

5.2. N — 3 with equal viscosities p^ — p^ = p^ 
This equal- viscosity case was first discussed by Stone fc Brenner] (1996). Putting this 



special case p {1) = /i (2) = /i (3) into (3.271 and solving it with MATLAB's Symbolic Math 



Toolbox (MuPAD), we obtain the same solution as equation (8) in Stone & Brenner 



(199G). 



5.3. Navier-Stokes and inviscid cases 

In Appendix [C] we validate the generalized form of our instability analysis against pre- 
vious results for inviscid and/or Navier-Stokes problems. For example, we find identical 



results to Chauhan et al. (2000) for the N = 3 case of a two- fluid compound jet sur- 



rounded by air (with negligible air density and viscosity). 

5.4. Comparison with numerical experiments 
Deng et al. ( 2011[ ) studied axisymmetric capillary instabilities of the concentric cylindrical 



shell problem (N = 3) for various viscosity contrasts by solving the full Navier-Stokes 
equation via finite-element methods. [The Stokes equation is a good approximation for 
their model, in which the Reynolds number is extremely low (Re « 10~ 10 ).] In particular, 
they input a fixed initial perturbation wavenumber k , evolve the axisymmetric equations, 
and fit the short-time behaviour to an exponential in order to obtain a growth rate. With 
their parameters i?W = 60 pm, = 120pm, 7W = = 0.6N/m, p^ = 10 5 Pa-s, 
and pS 1 ^ = p^ — rjpS 2 ' (77 = 10~ 4 , 10~ 3 , • • • , 10 3 ), we compute the maximum growth 
rate <7 max for each ratio r/ via the equation det(M^) = 0. For comparison, we also 
compute the growth rate cr max (fco) for their fixed k = 7.9 x 10 3 m _1 . [Because numerical 
noise and boundary artifacts in the simulations will excite unstable modes at k ^ fco, it 
is possible that cr max and not er max (fco) will dominate in the simulations even at short 
times if the former is much larger.] The inset of figure [2] plots the wavenumber A: max that 
results in the maximum growth rate versus the viscosity ratio rj (/i^ 1,3 ) / p^). In figureJSJ 
we see that the growth rate obtained by Deng et aT] ( |2011[ ) (blue circles) agrees well with 



the growth rate <r max (A:o) predicted by linear stability analysis (blue line) except at large 
viscosity contrasts (77 ^> 1 or r\ -C 1). These small discrepancies are due to the well-known 
numerical difficulties in accurately solving a problem with large discontinuities. 



6. N = 3 examples 



In this section, we study the three-fluid (N = 3) problem. Three or more concentric 



layers are increasingly common in novel fibre-drawing processes (Hart et a/.|2002 


Kuriki 


et al. 


2004 


Pone et al. 


2006 


Abouraddy et al. 2007 


Deng et al. 


2008 


). By exploring 



a couple of interesting limiting cases, in terms of shell viscosity and shell thickness, we 
reveal strong connections between the N = 3 case and the classic N = 2 problem. 
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a(k Q ) from numerical simulation 
a(k Q ) from linear stability analysis 
from linear stability analysis 




10 10 10 10 10 
viscosity ratio contrast: \r ' ' I n' 21 



Figure 2: Comparison between linear stability analysis and numerical experiments [data 
from Deng et al. (2011 1] for TV = 3 cylindrical-shell model. The growth rate <x max (fco) 



computed by Deng et al. (2011 ) numerically via the finite-element method (blue circles) 



agrees well with the growth rate predicted by linear stability analysis (dashed blue line) , 
except for small discrepancies in the regime of large viscosity contrast where accurate 
numerical simulation is difficult. The red line indicates the maximum growth rate cr max 
obtained by linear stability analysis. In the inset, the red line shows the wavenumbers 
fcmax for various viscosity-ratio contrasts and the dashed blue line represents the fixed 
fco used in numerical simulations. Model parameters: RS 1 ' = 60 fim, R (2 1 = 120^m, 



7 (1) = 7 ( 2 ) = 

and kn 



0.6N/m, = 10 5 
7.9 x ltfm- 1 . 



Pa-s, 



It 



(1) = ^(3) = 



r]fi 



(2) 



(rj = 10- 4 ,10" 



6.1. Case N = 3 and fi^ //p- 1 ' 3 ^ — > 0: shell viscosity <C cladding viscosity 
We first consider the limiting case in which the shell viscosity fj,^ is much smaller than 
the cladding vis cositi es fj,^ and fi^. Substituting (j,^ 2 '/^ 1 ' = and jjS 2 ^ / fi^ 3 ' 1 = into 
equation (3.271 gives simple formulae for the growth rates 

7 W (kR^) 2 -l 

ai{ ) "2m(W) 1 + ( ^( 1 ))2_ ( ^( 1))2 |^ ( • } 

and 

7< 2 > l-(fci?( 2 >) 2 

a2{ } ~ 2 ^ (3) * (2) 1 + (kR^Y ' 

Note that <7i(fc) is independent of 7^, R^ 2 \ and fi^ 3 \ while <J2{k) is independent of 7' 1 ), 
RS l \ and fi^. In particular, these growth rates are exactly the single-cylinder results 
predicted by Tomotika's model, as if the inner and outer layers were entirely decoupled. 
This result is not entirely obvious, because even if the shell viscosity can be neglected, it is 
still incompressible and hence might be thought to couple the inner and outer interfaces. 



Deng et al. (2011 1 conjectured a similar decoupling, but only in the form of a dimensional 



analysis.] 

6.1.1. Case N = 3, ^ 2 )/> (1 ' 3) -> and R^ -> 00 ^> N = 2 and /w/>in -> 

In the regime that the shell viscosity fi^ is much smaller than the cladding viscosities 
/iW and ^ 3 \ we further consider the limit i?( 2 ' — > 00. It corresponds to the case N = 2 
with a high viscous fluid embedded in another low-viscosity fluid, which must of course 
correspond exactly to Tomotika's case. From the asymptotic formulae of modified Bessel 
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(a) 

Figure 3: (a) Sketch of a very thin shell in 
R> X >[1 + e), surface-tension coefficients 7W 
£, we obtain an equivalent N = 2 geometry 

7 (l) +7 (2). 




a three-layer structure with radius R^ = 
and 7W. (b) In the limit of infinitesimal 
with a modified surface-tension coefficient 



functions Kq(z) and K 1(2) for large arguments ( Abramowitz fc Stegun|1992 ), we obtain 

7 ( 2 )|/fc| 



<7 2 (A) 



2/iO) 



< 



MS 



i? (2) -> +OO. 



(6.3) 



Thus, the growth rate of possible unstable modes is given by CTi(fc) in equation (6.1). 
Tomotika discussed this limiting case (N = 2) and gave a formula (37) ( |Tomotika|1935 ), 
which is exactly (6.1 1. 



6.1.2. Case N = 
The limit R*- 1 



' — > is equivalent to N = 



-s- <^=> N = 2 and /U ou t / Mm — > 00 
2 with a low-viscous fluid embedded in a 



high- viscosity fluid. For this case, it is easy to check that 02 (A; ) in equation (6.2) agrees 



with formula (36) in Tomotika (1935). However, we still have another unstable mode 



with a growth rate <Ji(k). Using the asymptotic formulae for Iq(z) and Ii(z) with small 
arguments ( Abramowitz &i Stegun|[l992 ) , we obtain 

^ (6 ' 4) 

This extra unstable mode o~\(k) results from the i nsta bility of a viscous cylinder with 
innnitesimally small radius R^\ In other words, (6.4) is the growth rate of a viscous 



cylinder in the air with a tiny but nonzero radius, which is also given by equation (35) 



of Lord Rayleigh (1892) 



that is, RW = 
This is moti- 



3 a 



6.2. Thin shell case: R (2 ^> = (1 + e), e -> 

Next, we study a three-layer structure with a very thin middle shell 
R^\\ + e) with e — > 0. A sketch of such a geometry is given in figure 
vated by a number of experimental drawn-fibre devices, which use very thin (sub-micron) 
layers in shells hundreds of microns in diameter in order to exploit optical interference 
effects ( |Hart et aL||2002| |Kuriki et a^|2004| |Pone et a/.||2006[ ) . 

Considering e as a small parameter, we expand the determinant equation ( 3.27[ ) in pow- 
ers of e. For a given wavenumber k, the two roots of this equation are a + (k) = ao(k)+0(e) 
and <r~(k) = 0(e 2 ), where ao(k) can be c ompu ted analytically by dropping the terms of 
order 0(e 2 ) in the determinant equation (3.27). After some algebraic manipulation, we 
find that cr (fc) actually is the solution for the N = 2 structure (i.e., ignoring the thin 
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relative thickness of shell layer e 
(a) In-phase mode a + 



relative thickness of shell layers 

(b) Out-of-phase mode a~ 



Figure 4: Two modes a + and a for thin-shell layer geometry, with = 1, R^> = 
RW(1 + e),jM = i, 7 (2) = 2, Ai ( 1 ) = 1,^( 2 ) = 2, = 3 and k = 0.5. (a) In-phase mode 
a + illustrates that the growth rate of the in-phase mode a + (k) for N — 3 approaches 
to the growth rate of N = 2 structure with the summed surface-tension coefficients as 
e — > 0. (b) Out-of-phase mode cr~ demonstrates that the out-of-phase growth rate cr~(k) 
decreases like e 2 as e — > 0. 



shell) with a modified surface-tension coefficient 7 



(i) 



v(2) 



It is also interesting to con- 



sider a limit in which grows as e shrinks. In this case, we find the same asymptotic 
results as long as (j,( 2 > /f/,^ 1 ' 3 ' grows more slowly than 1/e. Conversely, if it grows faster 
than 1/e, then the thin-shell fluid acts like a "hard wall" and all growth rates vanish. 
Instead of presenting a lengthy expression for er + (fc), we demonstrate a numerical ver- 
ification in figure |4] As indicated in figure 4a the growth rate a + {k = 0.5) for N = 3 



approaches the the growth rate of N 

e 0. The parameters are R^ 1 ' = = l,7 v ~' = = l, ju 

In figure 4b we show that the growth rate a~ (k) decreases like e 2 as e — > 0. 



2 with the summed surface-tension coefficients as 
l, 7 (i) = 1, 7 (2) = 2, M « = 1, M ( 2 ) = 2 and ^ = 3. 



To better understand these two modes, we consider the eigen-amplitudes at the two 
interfaces. For the mode with growth rate a + (k), two interfaces are moving exactly in 
phase. Since the thickness of this shell is so thin, it is not surprising that one can treat 
two interfaces as one with a modified surface-tension coefficient +7^) for this mode 
(see figure 3b and the inset of figure 4a). The eigen- amplitude (defined in section 3.5) 
corresponding to this in-phase mode is (1/2, 1/2), independent of 7^ and For the 
other mode, with growth rate a~(k), the two interfaces are moving out of phase (see 
the inset of figure [4b| . The eigen-amplitude for this out-of-phase mode is found to be 



(— 7^,7^)/(7^ + 7 ( - 2 ' ) ), which means that the two interfaces are moving in opposite 
directions with amplitudes inversely proportional to their surface tensions. Due to the 
tiny thickness of the shell compared to its radius of curvature, this case approaches the 
case of a flat sheet, which is known to be always stable (Drazin & Reid 2004), as can be 
proved via a surface-energy argument. 

A related thin-shell problem was investigated by Chauhan et al. ( 2000 ) for the Navier- 
Stokes equations with an inviscid (gaseous) outer fluid. Those authors also found that the 
problem reduced to N — 2 instabilities (single fluid surrounded by gas) with a summed 
surface tension. 
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Figure 5: Maximum effective growth rates versus wavenumber k. For a three-layer struc- 
ture with i?W = 1, i?( 2 ) = 5, [j,^' = a^ 2 -* = 1, and 7W = 1, the maximum effective growth 
rates tr^ x (A;) are plotted for several values of j( 2 K For — 12.19 (magenta line), there 
are two equal maximum effective growth rates cr^ x (fci ss 0.58) = cr°* x (fc2 ~ 0.114). 



7. Effective growth rate and competing modes 

In previous work on linear stability analysis, most authors identified the maximum a 



with the dominant breakup process (Lord Rayleigh 1879 1892 Tomotika 1935 1. This 



exclusive emphasis on the maximum a was continued in recent studies of N = 3 systems 
( |Stone fc BrennerjT9 96 ; Chauhan et al. 2000[ ), but here we argue that the breakup process 
is more complicated for TV > 2. In a multi-layer situation, however, there is a geometric 
factor that complicates this comparison: not only are there different growth rates a, but 
there are also different length scales over which breakup occurs. As a result, it is 
natural to instead compare a breakup time scale given by a distance (R) divided by 
a velocity, where R is some average radius for a given growth mode (weighted by the 
unit-norm eigen-amplitudes SR 1 -- 71 * 1 ). In our case, we find that a harmonic-mean radius R 
is convenient, and we define an effective growth rate (~ 1/breakup time ~ velocity/i?) 

by 



\k)=a. ] (k) 



n—1 



in) 



(7.1) 



Now, it is tempting to wonder what happens if two different wavcnumbers k\ and 
k,2 have the same maximum effective growth rates, a question that does not seem to 
have been considered in previous linear stability analyses. Let us consider a particular 



1,RW = 5, At 



(i) 



(2) 



1, and 7 



(i) 



1. The 



three-layer structure with R^ 

maximum effective growth rates o-^ ax (k) = maxj[cr| ff (/c)] versus k are plotted in figure [H] 
for several values of 7^. For example, at 7^ = 12.19, we find that of^^fci sa 0.58) = 
a max(^2 ~ 0.114), so that there are two competing modes at very different length scales 
2n/ki — 10.83 and 2n/k2 — 55.12. In contrast, for 7^ = 6 we see that the short length- 
scale instability should dominate, while at 7" = 25 the long length scale instability 
should dominate. 

To test our predictions, we implemented a full three-dimensional Stokes-flow numerical 
scheme to simulate the breakup process of this cylindrical-shell system. A brief description 
of this hybrid scheme, a combined spectral and level-set method, is given in Appendi x [B] 
We use initial white-noise perturbations on both interfaces R^ and R^ (see figure paE 
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(a) (b) 7 (2) = 6 (c) 7 (2) = 15 (d) 7 (2) = 25 



Figure 6: Numerical Stokes-flow simulations for three-layer systems with different 
(a) Initial white-noise perturbations of the interfaces. As predicted by maximum effective 
growth rates, the systems with 7^) = 6 (b) and 7 ( 2 - ) = 25 (d) exhibit breakup initially 
via the short- and long-scale modes, respectively (which are dominated by motion of the 
inner and outer cylinders, respectively). Near-simultaneous breakup occurs for 7W = 15 
(c). 



rather than u^ax would lead one to predict that 
4.15, in which case all three values in figure [6] 



The computational box is 16 x 16 x 108 with resolutions 160 x 160 x 480 pixels. As 
predicted, 7 ^ 2 ^ = 6 and 7 ' 2 ^ = 25 exhibit breakup initially via the short- and long-scale 
modes, respectively (which are dominated by motion of the inner and outer cylinders, 
respectively). It is interesting to estimate the intermediate 7 ^ 2 ^ where the two breakup 
processes occur simultaneously. Linear stability analysis predicts 7^ 2 -* ps 12.19, and indeed 
we find numerically that near-simultaneous breakup occurs for 7 < - 2 - ) ps 15 (see figure pq) . 
In contrast, simply looking at <r max 
simultaneous breakup occurs at 7W p 

would have looked like figure 6d (large scale dominating). In the case of near-simultaneous 
breakup timescales, the dominant breakup process may be strongly influenced by the 
initial conditions (i.e. the initial amplitudes of the modes), which offers the possibility of 
sensitive experimental tunability of the breakup process. 

The final breakup of the fluid neck is described by a self-similar scaling theory for the 
case of a single fluid jet (Eggers 1993), and so it is interesting to examine numerically 
to what extent a similar description is possible for the N — 3 system. In particular, at 
the last stage of a single-cylinder breakup process, a singularity develops at the point 
of breakup which does not possess a characteristic scale, and hence a set of self-similar 



profiles can be predicted ( Eggers & Villermaux 2008 ) . For both a viscous jet in gas ( Eggers 
1993 ) and a viscous thread in another viscous fluid (N = 2) ( jLister fc Stone|1998[1Cohen1 



et al. 1999| ), these principles predict that the neck radius h(t) vanishes linearly with 
time as h(t)/j,/j ~ (to — t) where to is the breakup time. However, there is no available 
scaling theory for TV > 3 systems. Here, we simply use our numerical simulations above 
to study the rate at which the neck radius vanishes in an N — 3 system. For all three 
cases with different 7^ 2 ', the neck radius of the outer interface vanishes with time in an 
asymptotically linear fashion as the breakup time is approached (see figure [7]). This is 
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Figure 7: Radius of the fluid neck versus time during the final phase of the breakup of the 
outermost interface, from the 3-dimensional Stokes simulations of figure [6] This breakup 
is asymptotically linear with time, similar to the predictions of the scaling theory for 
N = 2 systems ( jLister fc Stone|1998| [Cohen et aL||1999 |. 



not surprising in the 7^ 2 - 1 = 6 case where the inner surface has already broken up — the 
breakup of the outer surface reduces to an N — 2 problem when the neck becomes thin 
enough — and we find that h(t) \x / V 2 - 1 ~ 0.024(in — t), in reasonable agreement with the 
0.033 value predicted analytically ( |Cohen et aL|[l999| given the low spatial resolution 
with which we resolve the breakup singularity (h/RS^ =0.1 corresponds to 5 pixels). 
Moreover, we find that in this equal- viscosity N = 3 system, all three 7^) values yield 
slopes of h{t) jj, / that are within 10% of one another, indicating that the inner-surface 
tension = 1 has a relatively small impact on the outer-surface breakup. 



8. iV-Layer structures 

Since all previous work has studied only N — 2 or N — 3, it is interesting to consider 
the opposite limit of N — > 00. We consider two examples: a repeating structure of two 
alternating layers, and a structure with continuously varying viscosity, both of which are 
approached as N — > 00. In fact, concentric-shell structures with dozens of alternating 



fluid layers have been used experimentally in optical fibres (Hart et al. 2002 Kuriki 
et aL|2 004). However, the motivation of this section is primarily exploratory, rather than 
engineering — to begin to discover what new phenomena may arise for large N. 



First we consider an ABABAB ■ 



.1. Alternating structure 

structure of two alternating, repeating layers A and 



B as shown schematically in the left inset of figure [8j We choose \x 
jjy 1 } = 2 otherwise. The other parameters are = 1 + 0.2(n - 



(") = 1 if n is odd and 
1) and 7 (") = 1. 
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number of layers 



Figure 8: Growth rates of an ABABAB ■ ■ ■ alternating structure. Both er max (black line) 
and (red line) converge to finite asymptotic values as N — > oo, although in the 

former case the asymptotic value depends on the parity of N. The differences between 
""max ( re d line) and (<r max ) cff (blue line) imply that the modes corresponding to <r max and 
^max are n °t always the same. The right inset shows the eigen-amplitudes <5i? max (black 
dots) and <$-R^f ax (red circles) for N — 70, corresponding to er max and cr^ x respectively. 
The left inset depicts the structure whose parameters are = 1+0.2(71-1), 7 (n) = 1, 
and /i^™) = 1 if n is odd or fiS") = 2 otherwise. 



For this multilayer structure, we fi nd t hat both the maximum growth rate cx max of the 
fastest-growing mode and equation (7.1 1's maximum effective growth rate a^ax (corre- 
sponding to the shortest breakup time scale) apparently converge to finite asymptotic 
values as N — > oo (figure]!]). (We have checked that the absolute value of the slope of a max 
is monotonically decreasing for a broader range N values up to N = 120, and the slope 
is ~ 10~ 5 for N = 120. A rigorous proof of convergence requires a more difficult analysis, 
however.) The oscillations in figure]!] are due to the varying viscosity of the ambient fluid, 
which depends on the parity of N. It is interesting to know whether the fastest-growing 
mode and the mode with maximum effective growth are identical for a given N. To see 
this, we plot the effective growth rate of the fastest-growing mode (cr max ) cff versus N and 
compare it with cr^x versus N in figure [8j Note that (cr max ) cff = {max^fc[<jj(fc)]} cff and 
°max = max j,fe[ (7 | ff (&)]• From figurepl we can see that (er max ) cff and cr^ x are different for 
large N, which implies that the modes corresponding to cr max and <7^ x are not always 
the same. 

In the right inset of figure [ij we plot the eigen-amplitudes 8R max and (5ii°f ax for 
N = 70, corresponding to c max and cr^ ax respectively. The mode cr max is mostly motion 
of outer interfaces, while the mode <r^ x is mostly motion of inner interfaces. The mostly 
outer-interface motion for cr max explains why the value of cr max oscillates depending on 
the ambient fluid. Physically, the association of a^nx with the inner interfaces makes 
sense because, in our definition (7.1) of effective breakup rate, it is easier to break up 
at smaller radii (a smaller distance to breakup). Alternatively, if we defined "breakup 
distance" in terms of the thickness of individual layers, 
sense. 
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Figure 9: Growth rates for a continuous TV-layer structure. Both cr max and cr°f ax approach 
constants as N — > oo. The right inset plots the corresponding eigen-amplitudes <5i? max 
and SR!^ ax (for N=70). The left inset sketches the iV-layer structure: radius = i?i„, 



JV-l) 



(n) 



Rout, R^ n ' 



? (iV-l)_ 



N-2 



-(n — 1), and viscosity /Lt 



Min, yU 



(AT) 



/^out i 



and /Lt 
three-layer viscosity. 



AT-l 



(n — 1), approximating a continuously and linearly varying 



8.2. N -Layer structure for a continuous model 

In this subsection, we build an iV-layer model to approximate a three-layer structure 
with a continuous viscosity. The viscosity of intermediate layer /i m ;d of this three-layer 
structure is continuously varying from the viscosity of inner core /z; n to the viscosity of 
ambient fluid /i ut- A simple example is the linearly varying /^ m id, namely, /i m id('") = 
fi- m + ^"IZfi" i r — Rm), where i?j n < r < R ou t- The iV-layer structure (left inset of 
figure [9| with radius RW = i? inj R( N -V = R out , #(«) = fltt) + r(N ~^_- rW (n - 1), 

and viscosity fi^ — (j,i n , = /i ou t, and fi^ = fj,^ + M N Zi (n — 1) approaches 

this continuous model for large N. In order to obtain a physically realistic continuous- 
viscosity model with an energy that is both finite and extensive (proportional to volume) , 
we postulate a volume energy density ry CO nt analogous to surface energy. We approximate 
this by an iV-layer model constructed to have the same total interfacial energy: 



/•■Rout N ~ l 

/ 27rrdr = l [n) ^R {n \ (8.1) 

JRin „ =1 



?7cont 



where 77 cont is an appropriate energy (per unit volume) of the inhomogeneity. Correspond- 
ing to a uniform rj covA , the surface-tension coefficient in this N- layer structure is same 



on all the interfaces: namely 7W = j(N) for all n. From (8.1 1, we obtain the equivalent 
surface tension 7 (AT) in an iV-layer structure 
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With the parameters described above, we compute the maximum growth rate <7 max and 
the maximum effective growth rate <r^ x for this TV-layer structure. As shown in figure [9J 
both fj max and cr^ x approach constants as N — > 00, which should be the corresponding 
growth rates of the continuous three-layer model. In this example, o-^ ax and (cr max ) eff are 
the same for all N. The corresponding eigen- amplitudes 6R max and 5R^ ax (for iV=70) 
are plotted in the right inset of figure [9j 



9. Conclusions 

In this paper, we presented a complete linear stability analysis of concentric cylindrical 
shells in the Stokes regime (with the Navier-Stokes regime in AppendixjC]) and considered 
a few interesting examples and limiting cases. Many possibilities present themselves for 
future work. First, even in the cylindrical Stokes regime, only a few combinations of thick- 
nesses and material properties have been considered so far — it seems quite possible that 
consideration of larger parameter spaces, perhaps aided by computational optimization, 
could identify additional regimes for breakup processes, such as competitions between 
additional length scales or "effective" properties in many-layer systems that differ sub- 
stantially from the constituent materials. Second, one could extend this work beyond the 
incompressible Navier-Stokes regime to include fluid compressibility or even other phys- 
ical phenomena such as viscoelasticity that may play a role in experiments (for example, 
fibres are drawn under tension). Third, one could consider non-cylindrical geometries. 
This seems especially important in light of the recent experimental observations of az- 



imuthal breakup in cylindrical thin-shell fibre structures (Deng et al. 20081, since this 
paper points out that azimuthal breakup cannot arise in purely cylindrical structures 
(at least, not from surface tension alone). Instead, one may need to consider the "neck- 
down" structure of the fibre-drawing process, in which a large preform is pulled to a 
long strand with a much smaller diameter. More generally, such intriguing experimental 
results indicate that a rich variety of new instability phenomena may arise in emerging 
multi-fluid systems, with corresponding new opportunities for theoretical analysis. 
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Appendix A. Computations of the curvature 

Here we derive the curvature terms in (3.20). The level-set function (p n '{r, z,t) = 
corresponds to the n-th interface. The unit outward normal vector of this interface is 



(n 



(™) „(")'* 



v dr ' dz 



(1, -ifc&RWe^*-'"*)) 



|V0(")| ^/ ( O0W )2 | v/l + 0[(^i?(")) 2 ] 

and the unit tangential vector is 

(-ikSRW&( k *- ui \-T) 



= {n{ n \ -n^ = 



(Al) 



(A 2) 
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The curvature can now be computed as 

dn[ n) 



1 



dr r 



dz 
1 



/l + 0[pW)2] R(n) + ^WeK^-*) ^ + Q [{5R{n)? 



O 



{5R (n) f 



R{n) 



SR (n) 



(Rwy 



e i{kz- U t) + q W§ R (n)^ 



(A3) 



The normal velocity of the fluids on the interface must be equal to the normal velocity 

of that interface, and thus dC g t — vS n ^ ■ n n on the interface r — ^ n \z,t), where 

is the velocity vector. For the at-rest steady state (3.3) and ( |3.4[ ), to the lowest order in 

SR, this gives 

(A 4) 



\uj8R {n) ^5u^\R [n) ). 



Note that (A4) establishes the relation between the displa ceme nt amplitude 5R> n ' and 
the interface velocity Sui n \R^). Substituting (A4| into (A3), we obtain the lowest- 
order curvature in terms of the interface velocity Sui n \R^): 

1 



R(n) 



(i?(™)) 2 



^i(kz— tot) 



(A 5) 



Appendix B. Full 3-dimensional Stokes-flow numerical simulation 
scheme for coupled cylindrical-shell system 

In this section, we briefly present the numerical scheme that we used in §[7] to simulate 
the instabilities of coupled cylindrical-shell systems. We adopt a 3-dimensional Cartesian 
level-set approach. We use a separate level-set function cj)^ to denote each interface, and 
generalize the formulation of Chang et al. ( 1996 ) to N fluids by using N — 1 level-set 
functions governed by the following equations: 



N-l 



Vp + V- MV[/ + VC/ T )| =^ 7 (™) < 5(^))^("))-J-, (Bl) 

n— 1 ' 



and 



;)i ,U-V^ n > =0, (B2) 

where U is velocity, p is pressure, = denotes the interface between the n-th and 
(n + l)-th layers, 

7 («) i s 

the surface-tension coefficient of the n-th interface, <$(■) is a Dirac 
delta function, k is curvature, and /i is viscosity. 

The viscosity /x, now defined in the whole coupled system, is 

N-x 

fi(x) = + - V in) )H(4> {n) (x)), (B3) 

71=1 

where H( ) is the Heaviside step function. The curvature k(0W) can be computed directly 
by = V • |^0(n)| } since r^^jj is the unit outward normal vector of the n-th 

interface. 




Time 



Figure 10: The aspect ratio of a 2-dimensional elliptical blob versus time, obtained by 
different methods and implementations. For the system initially bounded by a; 2 /4 + 
y 2 = 1, with the elliptical blob viscosity /x; n = 1, the ambient fluid viscosity ^ out = 



and the surface tension 7 = 1, Buchak (2010) used the conformal mapping method via 



finite element implementation, obtaining the black evolution curve. The evolution curves 
(green, blue and red) given by our simulations converge to the black curve with the 
increasing resolutions and as the ambient viscosity ^S 2 ' goes to zero. 



Given the level-set function <fi( n \x,t) at time t, we first solve the steady Stokes equa- 
tions (Bl) to obtain the velocity U(x,t). With the known velocity U(x,t) at time t, 
the level-set function <j)( n \x, t + At) can be obtained by solving the convection equation 

In our implementation, the computation cell is a box with dimensions a x a x I in Carte- 
sian coordinates, with periodic boundary conditions. We choose a and I large enough such 
that the periodicity does not substantially affect the breakup process. We solve equa- 
tions (B 1 ) by a spectral method: we represent U and p by Fourier series (discrete Fourier 
transforms). For the constant fi case of §[7J ( |B 1| is diagonal in Fourier space and can 
be solved in one step by fast Fourier transforms (FFTs). More generally, for variable 
viscosity, we find that an iterative solver such as GMRES (generalized minimal residual 
method) or BiCGSTAB (biconjugate gradient stabilized method) (Barrett et al. 1994) 
converges in a few iterations with a constant /i preconditioner (i.e., block Jacobi) using 
the average [i. The level-set functions are described on the same grid, but using finite dif- 
ferences: WENO (weighted essentially non-oscillatory: Liu et al\ ( 1994[ )) in space, third 
order TVD (total variation diminishing) Runge-Kutta method in time ( |Shu fc Osher 



1989). The S(-) function is smoothed over 3 pixels with a raised-cosine shape (Osher & 



Fedkiw 2002). We use the reinitialization scheme of Sussman et al. (1994) to preserve 
the signed distance-function property |V<p"-*| = 1 of ^■ n 1 after each time step. 

Our simulation code is validated against a well-studied case: the evolution of a 2- 



dimcnsional elliptical blob (Kuiken 1990 Hopper 1991 Tanveer & Vasconcelos 1995 



Crowdy 2002 2003). It is known that the plane Stokes flow, initially bounded by a 



simple smooth closed elliptic curve, will eventually become circular under the effect of 
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surface tension. Crowdy (2002) illustrated that the evolution via a series of ellipse shapes 
is remarkably good approximation to the dynamics of a sintering ellipse [even though 
Hopper ( 1991 ) showed that the exact evolution shapes are not strictly elliptical]. Suppose 
the plane Stokes flow is bounded by the ellipse x 2 /4 + y 2 = 1 at t — 0. The viscosity 
of the elliptical blob p,i n = 1, the viscosity of ambient fluid ^ out = 0, and the surface 



tension 7=1. Buchak (20101 implemented the conformal mapping method (Tanveer & 



Vasconcelos 1995 Crowdy 2002 2003) and computed the evolution of the boundaries. 



The aspect ratio (the major axis over the minor axis) of the ellipses versus time is 
plotted (black curve) in figure 10 Since our simulation code only works for non-zero 
p^ n \ the evolution under py-> = 1 and p^ — ^ is expected to converge to the black 
evolution curve obtained by the conformal mapping. In figure |10[ we also plotted the 
evolution curves from our method with p^ =0.1 and resolution 256 x 256 (green curve), 
p^ = 0.01 and resolution 256 x 256 (blue curve), and p( 2 > = 0.01 and resolution 512x512 
(red curve). With high resolutions and small p^ 2 \ the evolution curves obtained by our 
simulation codes converge to the one given by a different method with an independent 
implementation. 



Appendix C. 



Linear stability analysis for concentric fluid shells 
governed by the full Navier Stokes equations 



In this section, we extended our linear stability analysis to concentric cylindrical fluid 
shells governed by the full Navier-Stokes equations. Let p( n > and p^ denote the density 
and viscosity of the n-th fluid; uf- n (z) is the radial component of the velocity and u z (z) 
is the axial component of the velocity. Following a linear stability analysis similar to §|3.2[ 
we find that the pressure (z) in the n-th fluid still satisfies Laplace's equation (|3.7[). 



Therefore, the perturbed pressure still satisfies the modified Bessel equation (3.8) and 



the solution in (3.9) is still valid. The velocity is obtained by solving the linearized 



Navier-Stokes equations 
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du 
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(CI) 
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dt 



An) 
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2„,(") 1 



d u 



(C2) 



Note that the nonlinear convection terms do not appear in the linearized equations 



(CI) and (C2) because the basic steady state (3.3)-(3.4) is at rest. Substituting the 



of the perturbed velocity in (3.10) is replaced by 



perturbed pressure (3.9) into equations (C1)-(C2), we find that the radial component 



Su^ir) = 



{n )K x {kr) - ifi(kWr) {n) h(kr)-h{k^r) 



jp( n )/k 



3 2 M (»)fc 



(n) /l(fc ( "M 

^ 2 M (»)fc 



(C3) 
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and the axial component of the perturbed velocity in (3.11 1 now becomes 

, (n)( {n) K (kr)-X^K (k^r) (n) Ip(kr) - A<"> I (k™r) 

Uz [r) ~ Cl pM u /k 2 pWw/k 

[n) i\^K Q (k^r) ( »)iAWj (fc(")r) 
Ca 2/iW* 4 2/x(»)fc : 



where 



A" 



-iwjo(") 



/' 



(«)fc 2 



and 



k (n) = X (n) k _ 



(C4) 



(C5) 



After matching boundary conditions (3.14), (3.16), (3.19) and ( 3.21 ), we can obtain the 
dispersion relat ion b y solv ing th e same determinant equation (3.27), except that A^ n ' n ) 
and fiW from (3.23) and (3.24) are replaced by 



,(ra,n') 

1 ns 



p("')/2fc 2 



-Il[*< n '«(">] 
p("')/2fc 2 
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ii[fc< 7 ' / )fl<")] 
M< "' ) , 

^) 

A 4 



where 
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and 



d(«) _ 
D ns — 
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(fcij(")) 2 



p(")/2/c 2 







• ( C12 ) 

7 '")_j l[fc (") fl (")] gjJfcWfiW] j^fct")^")] 

p(")/2fe 2 p» p» 

Note that, because of the u in A 1 -™' and the matrix becomes nonlinear in uj 

(or l/ui), and can no longer be reduced to a generalized eigcnproblem. Instead, one must 
solve the nonlinear eigenproblem (&> w )£ = iw/W^ (fc, Numerous methods have 
been developed for such problems ( Andrew et «Z.|1995] Guillaume|199 9; Ruhc 2006 Voss 
2007||Liao et aL||2010[ ). 
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the limit p 



(») 



The formulae (3.23) and (3.24) of A^ n ' n > and fiW for Stokes flow can be obtained 



0. As mentioned in 



directly from the formulae ( C 6 1 and (C 12 1 of A^ n ^ and B^) for general flow by taking 



5.1 



the most straightforward formulation of the 
Navier-Stokes matrices yields dependent columns when p — > 0. Here, we have chosen an 
appropriate linear combination of columns to avoid this difficulty, which is equivalent to 



the procedure suggested by Tomotika ( 1935 1 



However, the corresponding formulae for inviscid fluids cannot be obtained by simply 
taking the limit — > 0. For any small but nonzero //™\ the current formulation 



takes into account the boundary-layer effects (Batchelor 1973) by imposing a no-slip 



condition. For inviscid flows, we cannot assume that the axial velocities are continuous 
across the interfaces, since no-slip boundary conditions are not applied. Whenever the 
no-slip boundary condition is not applied, one has an additional degree of freedom, the 
equilibrium-state velocities Uz of the layers. This is easily incorporated, because it 

-in) 

merely converts several oj expressions to co — u z k. This happens in two places. First , 
the velocity adds an additional inertial term -p^u^dul^/dz to the left side of (|Cl) 



and -pWui"'^- to the left side of (|C2|). Second, it adds a new first-order term to 



equation ( 3.13 ) for continuity of no rmal velocity, since there is a term from multiplied 
by the 5R [n > in the numerator of (|A l|) for the n orma l vector. These terms change u) to 

z™ in (C7) for oS- n \ and they also multiply 
i^k) [cancelling the 1/lu factor multiplying 



k in (C5) for and to 



every B ^J mat rix (including 6 n s ) by lu/(lu 



M. 



(JV) 



in (3.26) 



the A^ n ' and matrices are otherwise unchanged, since u 



(n) 



must 



be equal for adjacent viscous layers. If the n-th fluid is inviscid while the (n— l)-th and/or 



(n + l)-th fluids are viscous, then /4 ns 



(n-l,n-l) 



and 6 n g 1 \ and/or A\^g 



(n,n+l) 



respectively, 



become 3x4 matrices that can be obtained from ( C 6 ) and \C 12 ) by eliminating the 
second row (corresponding to continuity of the tangential component of the velocity). If 

the n-th layer is inviscid, regardless of the adjacent layers, A. ' n \ and B^ , are 3x2 

•> . . . mvsd invsq 

matrices: not only has continuity of the tangential component of the velocity disappeared, 



but also the d/dr derivatives of the velocities in the momentum equations (CI) and ( C 2 1 



disappear when p^ — 0, eliminating the 4 j degrees of freedom. (This eliminates the 
need for the linear combinations of columns mentioned above, further simplifying these 



matrices.) More explicitly, the A^ n ' n \ and B. w , matrices for an inviscid n-th layer are 
' mvsd invs d 



(») 



obtained by matching the boundary conditions (3.13), (3.19), and (3.21), giving 
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For the special case uf^ = (at- rest steady state), the above formulae are equivalent 
to the ones obtained by eliminating the second row, eliminating the thir d and fourt h 
columns, then taking the limit pS n > — > in the general formulae (C6) and (C12). 



More precisely, for the viscosity term in (CI) and ( C 2 ) to be negligible, one must have 



py 1 ! up( n '/k . (Note that this is length- and time-scale dependent, so the validity 
of neglecting viscosity terms depends on the uj and k of the dominant growth mode.) It 
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is also interesting to consider a Galilean transformation in which a constant v z is added 
to uf^ for all n, which cannot change the physical results. Here, because all uj factors 
are accompanied by — u z k, it is clear that such a transformation merely shifts all of 
the mode frequencies ujj(k) by v z k, which does not change the growth rates (the imagi- 
nary part), while the shift in the real frequency is simply due to the frequency-fc spatial 
oscillations moving past any fixed z at velocity v z . 

As a validation check, we find that our formulation gives the same dispersion relations 
for various Navier-Stokes cases discussed in the previous literature: e.g., a single inviscid 
jet in air (ignoring the air density and viscosity) (Lord Rayleigh 1879), a single viscous 
jet in air (ignoring the air density and viscosity) (Lord Rayleigh 1892 1, a single viscous jet 



with high velocity in air (considering the air density but ignoring the viscosity) (Sterling 



& Sleicher 


1975| and 


(Chauhan et al. 


2000) 



and a compound jet in air (ignoring the air density and viscosity) 
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